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ABSTRACT 

This  Final  Report  describes  our  program  in  studies  of  (laminar/turbulent)  stability  and 
transition  in  non-equilibrium-chemistry  flows  characteristic  of  those  on  the  forebodies  of 
hypersonic  vehicles.  The  configuration  best  modelling  a  hypersonic  vehicle  is  an  elliptic  cone. 
Specifically,  we  investigated  and  optimized  a  Parabolized  Navier-Stokes  solution  for  the  basic-state 
flow  past  a  sharp  elliptic  cone  including  the  region  between  the  wall  and  the  shock.  We  formulated 
the  Parabolized  Stability  Equations  for  3-D  flows. 
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1.  Introduction 

This  Final  Report  describes  our  program  in  studies  of  (laminar/turbulent)  stability  and 
transition  in  non-equilibrium-chemistry,  bounded  flows  characteristic  of  those  on  the  forebodies  of 
hypersonic  vehicles.  In  this  Report,  Section  2  contains  a  list  of  accomplishments.  Section  3 
contains  a  summary  of  the  work. 

2.  Technical  Accomplishments 

In  the  past  3  years,  2  students  have  been  supervised,  15  publications  have  been  written  or  are 
in  preparation,  and  9  talks  and  lectures  have  been  given. 


Publications 

"Transition  Correlations  in  3-D  Boundaiy  Layers,"  H.L.  Reed  and  T.S.  Haynes,  AlAA  Journal, 
Volume  32,  Page  923,  1994. 

"Linear  Disturbances  in  Hypersonic,  Chemically  Reacting  Shock  Layers,"  G.K.  Stuckert  and  H.L. 
Reed,  AIAA  Journal,  Volume  32,  Number  7,  Page  1384,  1994. 

"Lineal'  Stability  Theory  Applied  to  Boundary  Layers,"  H.L.  Reed,  W.S.  Saric,  and  D.  Arnal, 
Annual  Review  of  Fluid  Mechanics,  Volume  28,  1996. 

“CFD  Validation  Issues  in  Transition  Modelling,”  H.L.  Reed,  T.S.  Haynes,  and  W.S.  Saric, 
submitted  to  AIAA  Journal. 

“Numerical  Investigation  Nonlinear  Saturation  of  Crossflow-Dominated  Boundary  Layers  Using 
PSE,”  T.S.  Haynes  and  H.L.  Reed,  submitted  to  Journal  of  Fluid  Mechanics. 

“Investigation  of  PSE  Initial  Conditions  for  Crossflow-Dominated  Boundary  Layers,”  T.S. 
Haynes  and  H.L.  Reed,  to  be  submitted  to  Journal  of  Fluid  Mechanics. 

“Spatial  Direct  Numerical  Simulations,”  H.L.  Reed,  to  be  submitted  to  Progress  in  Aerospace 
Sciences. 

"Dkect  Numerical  Simulation  of  Transition;  The  Spatial  Approach,"  H.L.  Reed,  Invited  Paper, 
AGARD  Course  in  Transition  Prediction  and  Modelling,  AGARD  Report  No.  793,  VonKarman 
Institute  and  Madrid,  March  1994. 

"Use  of  Transition  Correlations  for  Three-Dimensional  Boundaiy  Layers  within  Hypersonic 
Flows,"  I.J.  Lyttle  and  H.L.  Reed,  AIAA-95-2293,  26th  AIAA  Fluid  Dynamics,  Plasma 
Dynamics,  and  Lasers  Conference,  San  Diego,  June  19-22,  1995. 

"Use  of  Transition  Correlations  for  Three-Dimensional  Boundary  Layers  within  Hypersonic, 
Vi.scoLis  Flows,"  I.J.  Lyttle  and  H.L.  Reed,  Second  Symposium  on  Transitional  and  Turbulent 
Compressible  Flows,  1995  Joint  ASME/JSME  Fluids  Engineering  Conference,  Hilton  Head, 
August  1995. 

"Computations  in  Nonlinear  Saturation  of  Stationary  Crossflow  Vortices  in  a  Swept-Wing 
Boundary  Layer,"  T.S.  Haynes  and  H.L.  Reed,  AIAA  34th  Aerospace  Sciences  Meeting  and 
Exhibit,  AIAA-96-0182,  January  1996. 

"CFD  Validation  Issues  in  Transition  Modelling,"  T.S.  Haynes,  H.L.  Reed,  and  W.S.  Saric, 
Invited  Paper,  AIAA-96-2051 ,  Special  Session  on  Verification  and  Validation:  Uncertainties  and 
Special  Panel  on  Quantification  of  CFD  Uncertainties,  27th  AIAA  Fluid  Dynamics  Conference, 
New  Orleans,  June  17-21,  1996. 
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“Drag  Prediction  and  Transition  in  Hypersonic  Flow,”  H.L.  Reed,  W.S.  Saric,  R.  Kimmel,  S. 
Schneider,  and  D.  Arnal,  Invited  Paper,  AGARD  Inteipanel  (FDP&PEP)  Symposium  on  "Future 
Aerospace  Technology  in  Service  to  the  Alliance",  Paris,  France,  April  14-18,  1997 

“Drag  Prediction  and  Transition  in  Hypersonic  Flow,”  H.L.  Reed,  R.  Kimmel,  S.  Schneider,  and 
D.  Arnal,  Invited  Paper,  AIAA-97-I818,  Snowmass,  June  1997. 

“Role  of  Direct  Numerical  Simulations  in  Transition  Modelling,”  H.L.  Reed,  Invited  Paper,  First 
AFOSR  DNS/LES  Conference,  New  Orleans,  August  1997. 

Presentations 

"Direct  Numerical  Simulation  of  Transition;  The  Spatial  Approach,"  H.L.  Reed,  Invited  Paper, 
AGARD  Course  in  Transition  Prediction  and  Modelling,  AGARD  Report  No.  793,  VonKarman 
Institute  and  Madrid,  March  1994. 

"Use  of  Transition  Correlations  for  High-Speed  Three-Dimensional  Boundary  Layers,"  I.J.  Lyttle 
and  H.L.  Reed,  Bulletin  of  the  American  Physical  Society,  Volume  39,  Number  9,  Page  1913, 
November  1994. 

"Use  of  Transition  Correlations  for  Three-Dimensional  Boundary  Layers  within  Hypersonic 
Flows,"  I.J.  Lyttle  and  H.L.  Reed,  AIAA-95-2293,  26th  AIAA  Fluid  Dynamics,  Plasma 
Dynamics,  and  Lasers  Conference,  San  Diego,  June  19-22,  1995. 

"Use  of  Transition  Correlations  for  Three-Dimensional  Boundary  Layers  within  Hypersonic, 
Viscous  Flows,"  I.J.  Lyttle  and  H.L.  Reed,  Second  Symposium  on  Transitional  and  Turbulent 
Compressible  Flows,  1995  Joint  ASME/JSME  Eluids  Engineering  Conference,  Hilton  Head, 
August  1995. 

"Computations  in  Nonlinear  Saturation  of  Stationary  Crossflow  Vortices  in  a  Swept-Wing 
Boundary  Layer,"  T.S.  Haynes  and  H.L.  Reed,  AIAA  34th  Aerospace  Sciences  Meeting  and 
Exhibit,  AlAA-96-0182,  January  1996. 

"CFD  Validation  and  Verification  Questions  in  Transition  Modelling,"  T.S.  Haynes,  H.L.  Reed, 
and  W.S.  Saric,  Invited  Paper,  AIAA-96-2051 ,  Special  Session  on  Verification  and  Validation: 
Uncertainties  and  Special  Panel  on  Quantification  of  CFD  Uncertainties,  27th  AIAA  Fluid 
Dynamics  Conference,  New  Orleans,  June  17-21,  1996. 

“Drag  Prediction  and  Transition  in  Hypersonic  Flow,”  H.L.  Reed,  W.S.  Saric,  R.  Kimmel,  S. 
Schneider,  and  D.  Arnal,  Invited  Paper,  AGARD  Inteipanel  (FDP&PEP)  Symposium  on  "Future 
Aerospace  Technology  in  Service  to  the  Alliance",  Paris,  France,  April  14-18,  1997 

“Drag  Prediction  and  Transition  in  Hypersonic  Flow,”  H.L.  Reed,  R.  Kimmel,  S.  Schneider,  and 
D.  Arnal,  Invited  Paper,  AIAA-97-1818,  Snowmass,  June  1997. 

“Role  of  Direct  Numerical  Simulations  in  Transition  Modelling,”  H.L.  Reed,  Invited  Paper,  First 
AFOSR  DNS/LES  Conference,  New  Orleans,  August  1997. 


Ph.D.  Students 

T.  Haynes,  "Nonlinear  Stability  and  Saturation  of  Crossflow  Vortices  in  Swept-Wing  Boundary 
Layers,"  completed  Fall  1996. 

I.  Lyttle,  "Stability  of  Hypersonic  Flow  over  an  Elliptic  Cone,"  expected  Summer  1998. 


Reed:  Hypersonics  (AFOSR) 


page  6 


3.  Completed  Work 

3.1  Introduction 

The  world  over,  theorists  and  experimentalists  are  investigating  the  problems  associated  with 
hypersonic  flight.  Much  of  the  current  research  into  hypersonic  flight  is  concentrated  in  the  range 
of  Mach  5  to  Mach  15  (Reshotko,  1992).  The  skin-friction  drag  and  heat-transfer  rates  of 
hypersonic  vehicles  depend  on  the  state  of  the  boundary  layer,  that  is,  whether  it  is  laminai'  or 
turbulent.  Moreover,  the  perfoiTnance  of  the  engine  depends  on  the  state  of  the  boundary  layer  at 
the  inlet.  But  the  characteristics  of  transition  at  these  speeds  are  not  well  understood.  Transition  is 
known,  however,  to  be  highly  dependent  on  the  details  of  the  flowfield. 

The  current  conceptional  design  of  the  hypersonic  vehicle  includes  a  forebody  which  could  be 
modeled  as  a  shaip  cone  with  an  elliptical  cross-section.  It  is  this  geometry,  at  zero  angle-of- 
attack,  that  is  of  interest  to  this  study.  The  boundary  layers  associated  with  this  geometry  ai'e 
three-dimensional.  As  Reed  &  Saiic  (1989)  point  out,  three-dimensional  boundary  layers  ai'e 
susceptible  to  crossflow  instabilities,  as  well  as  streamwise  instabilities.  In  fact,  these  crossflow 
instabilities  are  often  the  dominant  mechanisms  responsible  for  transition,  and,  in  discussions  on 
upper-atmosphere  hypersonic  transition  experiments,  much  consideration  has  been  given  to  the 
minimization  of  three-dimensional  effects  (Kimmel,  1994). 

Hypersonic  flows  are  more  complicated  than  subsonic  and  supersonic  flows  for  some  of  the 
following  reasons.  1)  At  hypersonic  speeds,  the  gas  often  cannot  be  modeled  as  perfect  because 
the  molecular  species  begin  to  dissociate  due  to  aerodynamic  heating.  In  fact,  sometimes  there  ai'e 
not  enough  intermolecular  collisions  to  support  local  chemical  equilibrium  and  a  nonequilibrium- 
chemistry  model  must  be  used.  2)  The  bow  shock  is  very  close  to  the  edge  of  the  boundary  layer 
and  must  be  included  in  the  calculations.  3)  The  boundary  layers  on  the  forebody  are  highly  3-D. 
All  of  these  effects  must  be  included  in  predictions  of  transition. 

3.2  Chemically  Reacting  Basic-State  Flow 

The  purpose  of  this  investigation  is  to  examine  the  thi'ee-dimensional  nature  of  boundary 
layers  found  on  sharp  cones  of  elliptical  cross  section.  Initially,  the  basic  state  is  solved  using  a 
finite-difference  formulation  of  a  set  of  Parabolized  Navier-Stokes  (PNS)  equations.  The  set  of 
PNS  equations  used  is  that  derived  by  Lubard  &  Helliwell  (1974).  The  finite-difference  algorithm 
used  to  solve  these  PNS  equations  descends  from  the  scheme  derived  by  Tannehill  et  al.  (1982). 
Since  the  boundary  layer  occupies  a  substantial  fraction  of  the  shock  layer  at  hypersonic  speeds, 
these  equations  are  transformed  into  a  coordinate  system  where  the  basic-state  bow  shock  is  a 
boundary  of  the  computational  domain.  We  recognize  that  for  complex  flowfields  with  strong 
viscous/inviscid  interaction,  reduced  forms  of  the  equations  of  motion  do  not  provide  solutions 
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accurate  enough  for  transition  studies,  and  it  is  necessary  as  a  next  step  to  solve  either  the  full 
Navier-Stokes  (NS)  or  Thin-Layer-Navier-Stokes  (TLNS)  equations.  We  are  in  the  process  of 
formulating  and  completing  this  activity. 

Code  Optimization 

We  have  already  made  major  improvements  to  increase  the  computational  efficiency  of  the 
PNS  scheme.  The  reason  we  did  this  was  in  preparation  for  the  addition  of  the  nonequilibrium- 
chemistry  terms  and  the  extension  to  either  the  full  Navier-Stokes  or  Thin-Layer-Navier-Stokes 
formulations,  which  (we  knew  from  our  previous  experience)  would  significantly  increase  the 
magnitude  of  the  problem. 

For  the  Cray  Y-MP  and  the  Cray  C-90,  there  is  a  facility  called  vectorization,  whereby  the 
CPU  can  act  like  an  assembly  line.  This  allows  many  operations  to  be  performed  at  one  time, 
increasing  the  number  of  arithmetic  operations  performed  without  increasing  the  amount  of  CPU 
time  used.  To  achieve  vectorization,  certain  rules  must  be  followed  concerning  the  order  of 
operations,  etc.  In  many  cases,  this  requires  major  changes  in  the  algorithm  and  corresponding 
changes  in  the  code.  For  this  research,  vectorization  is  possible  and  is  exploited.  The  most 
computationally  intensive  parts  of  the  program  involve  the  solution  of  block-tridiagonal  systems  of 
equations.  The  majority  of  the  optimization  effort  was  placed  here. 

Step  1.  A  routine  was  written  to  invert  mati'ices  on  the  diagonal  to  exploit  their  structure. 

Step  2.  Major  input  and  output  routines  use  binary  files  instead  of  ASCII  files. 

Step  3.  The  block-tridiagonal  systems  are  solved  plane-by-plane  instead  of  line-by-line. 

Step  4.  The  algorithm  was  analyzed  to  remove  code  which  was  not  necessary. 

A  simple  representative  run  of  the  code  was  performed  to  gauge  the  effectiveness  of  the 
optimization.  The  speed  improvements  do  not  correspond  exactly  to  the  improvements  in  time 
because  Step  4.  was  carried  out  over  the  entire  course  of  the  optimization.  The  top  speed  of  one 
CPU  of  the  Cray  Y-MP  is  338  MFLOPS.  According  to  the  staff  at  CEWES,  this  is  a  theoretical 
limit  that  is  rarely  approached.  However,  they  are  supportive  of  efforts  to  achieve  such  goals. 


Code  Speed  (MFLOPS)  CPU  Time  to  complete  (s) 

Before  optimization  10.4  207 

After  Step  I  40.9  33 

After  Step  2  60.7  22 

After  Step  3  163.6  8 
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SIGNIFICANCE  One  obvious  benefit  of  this  effort  is  the  major  savings  of  CPU  resources 
for  this  project.  Further  more  complex  and  resource  intensive  codes  written  in  support  of  this 
research  will  be  written  using  the  principles  learned  in  this  simple  initial  effort.  This  will  result  in 
more  efficient  use  of  a  valuable,  limited  resource  to  solve  a  critical  national  problem. 

Chemistry 

Qualitative  as  well  as  quantitative  differences  in  the  transition  of  the  shock  layer  on  a  forebody 
can  be  observed  when  the  effects  of  finite-rate  chemical  reactions  are  included.  We  ai'c  following 
the  studies  of  stability  on  a  cone  of  circular  cross-section  by  Stuckert  &  Reed  (1994)  who  used  a 
five-component  model  for  dissociated  air:  O2,  N2,  NO,  N,  and  O.  Only  the  effects  of  dissociation 

were  considered;  those  of  ionization  were  not.  The  mixture  was  also  assumed  to  be  one  of  ideal 
gases  in  thermal  equilibrium.  The  viscosity  used  to  determine  the  viscous  stress  tensor  was 
computed  using  the  mixture  laws  of  Brokaw  (1958).  The  translational  thermal  conductivity  was 
computed  similarly,  whereas  the  internal  thermal  conductivity  was  determined  as  described  in 
Hirschfelder  (1957).  The  molar  fluxes  were  computed  using  the  multicomponent  diffusion  model 
described  by  Curtiss  &  Hirschfelder  (1949),  but  only  diffusion  due  to  concentration  gradients  was 
included  -  diffusion  due  to  pressure  and  temperature  gradients  and  body  forces  was  not.  Finally, 
the  law  of  mass  action  was  used  to  compute  the  molar  rate  of  production  of  each  species  assuming 
that  they  participate  in  the  elementary  reactions: 

02  +  M<->0  +  0  +  M 
N2-1-M0N  +  N  +  M 
NO  +  M<->N  +  0  +  M 
N2  +  O  NO  -I-  N 
O  -I-  NO  N  -I-  O2 

where  M  is  a  collision  pai'tner  (any  of  the  species  present  in  the  mixture)  which  transfers  energy  in 
a  reaction. 

The  thermophysical  data  needed  for  the  analysis  were  taken  from  a  variety  of  sources. 
Collision  cross  section  data  for  the  transport  properties  were  found  in  Biolsi  &  Biolsi  (1983), 
Biolsi  (1988),  Capitelli  &  Devoto  (1973),  Capitclli  &  Ficocelli  (1972),  Cubley  &  Mason  (1975), 
Levin  et  al.  (1987,  1988),  Monchick  (1959),  and  Yun  &  Mason  (1962).  Thei-modynamic  data 
were  taken  from  Blottner  et  al.  (1971)  and  Jaffe  (1987).  Reaction-rate  data  were  found  in  Camac 
&  Vaughan  (1961),  Wray  (1962),  Thielen  &  Roth  (1986),  Monat  et  al.  (1978),  and  Hanson  & 
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Salimian  (1984).  Reverse  reaction-rate  constants  were  computed  using  the  law  of  detailed  balance 
to  express  them  in  terms  of  the  forward-rate  and  equilibrium  constants.  For  a  complete  description 
of  the  constitutive  equations  as  well  as  detailed  references  for  the  thermophysical  data,  refer  to 
Stuckert  (1991). 

Josyula  &  Shang  (1993)  at  Wright  Patterson  Air  Force  Base  computed  the  2-D  steady,  basic- 
state  flow  over  a  hemisphere-cylinder  at  Mach  10-18  Hemisphere-Cylinder.  They  used  an  explicit 
NS  formulation  with  a  5-component  chemical  non-equilibrium  model,  comparable  to  Stuckert  & 
Reed's.  However,  the  objective  of  their  work  was  to  determine  the  sensitivity  of  heat- transfer 
predictions  to  the  temperature  model  used  (single  temperature  to  multi-temperature,  Stuckert  & 
Reed  used  single  temperature).  They  used  a  40  x  50  grid  on  a  CRAY  XMP,  which  performed  at 
3.38E-4  cpu-s/(pt. iteration)  and  required  20,000  iterations  to  converge. 

At  this  point  we  have  succcessfully  included  equilibrium  chemistiy  in  the  PNS  equations. 
Over  the  next  year,  we  shall  incoiporate  non-equilibrium  chemistry  and  investigate  the  effects  of 
(sensitivity  to)  chemistry  (equilibrium  and  non-equilibrium),  as  well  as  formulate  and  solve  either 
the  full  Navier-Stokes  (NS)  or  Thin-Layer-Navier-Stokes  (TLNS)  equations. 

3.3  Parabolized  Stability  Calculations  for  Transition  Prediction 
3.3.1  Introduction 

In  wall-bounded  shear  layers,  transition  occurs  because  of  an  incipient  instability  of  the  basic 
flow  field,  which  depends  intimately  on  subtle,  and  sometimes  obscure,  details  of  the  flowfield.  In 
other  words,  the  wall-bounded  shear  layer  is  an  open  system.  Disturbances  in  the  freestream,  such 
as  sound  or  vorticity,  enter  the  boundary  layer  as  steady  and/or  unsteady  fluctuations  of  the  basic 
state.  This  part  of  the  process,  called  receptivity  (Morkovin  1969,  Reshotko  1984),  provides  the 
vital  initial  conditions  of  amplitude,  frequency,  and  phase  for  the  breakdown  of  lamina:'  flow.  The 
recent  progress  in  this  area  is  summarized  in  Goldstein  &  Hultgren  (1989)  and  Saric  et  al.  (1994). 

Initially  these  disturbances  may  be  loo  small  to  measure  and  they  are  observed  only  after  the 
onset  of  an  instability.  The  initial  growth  of  these  disturbances  is  described  by  lineai'  stability 
theory  (LST).  Reed  et  al.  (1996)  review  this  subject. 

For  two-dimensional  (2-D)  boundary  layers,  this  growth  is  weak,  occurs  over  a  viscous 
length  scale,  and  can  be  modulated  by  pressure  gradients,  mass  flow,  temperature  gradients,  etc. 
As  the  amplitude  grows,  thiee-dimensional  (3-D)  and  nonlinear'  interactions  occur  in  the  form  of 
secondary  instabilities  (Herbei't  1988,  Cowley  &  Wu  1994,  Healy  1995).  At  this  point, 
disturbance  growth  is  very  rapid  (now  over  a  convective  length  scale)  and  breakdown  to  turbulence 
occurs  quickly. 
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3.3.2  Formulation 


In  this  Section  we  present  the  formulation  for  the  research.  Our  main  focus  is  the  stability  and 
transition  of  boundary-layer  type  flows.  For  transition  analysis  using  the  Parabolized  Stability 
Equation  (PSE)  approach,  equations  governing  the  disturbance  are  typically  solved  separately  from 
the  basic  state.  The  validity  of  the  basic-state  formulations  must  also  be  considered  since  the 
transition  process  is  known  to  be  sensitive  to  subtle  changes  in  the  basic  state.  In  all  cases  the 
numerical  accuracy  of  the  basic  state  must  be  very  high,  because  the  stability  and  transition  results 
will  be  very  sensitive  to  small  departures  of  the  mean  flow  from  its  “exact”  shape.  The  stability  of 
the  flow  can  depend  on  small  variations  of  the  boundary  conditions  for  the  basic  state,  such  as 
freestream  velocity  or  wall  temperature.  Therefore,  basic-state  boundai'y  conditions  must  also  be 
vei-y  accurate.  See  the  discussion  and  examples  of  Arnal  (1994)  and  Malik  (1990). 

3.3.3  Parabolized  Stability  Equations 


In  recent  years  the  parabolized  stability  equations  (PSE)  have  become  a  popular  approach  to 
stability  analysis  owing  to  their  elegant  inclusion  of  the  nonpaiallel  and  nonlinear  effects  which  ai-e 
ignored  by  LST  (Herbert  1994;  Bertolotti  1990).  For  linear  PSE  (LPSE),  we  consider  a  single 
monochromatic  wave  as  the  disturbance.  The  distui'bance  is  decomposed  into  a  rapidly  varying 
“wave  function”  and  a  slowly  vaiying  “shape  function”.  We  accomplish  this  here  with  a  multiple 
scales  approach. 

(p'{x,y,z,t)  =  ^{x,y)  ;K(x,z,r) +C.C.  (1) 

shcipejun  cli  on  wavefun  cli on 


where 


dx  ^  ^  dz  dt 


This  gives  the  following  form  for  the  streamwise  derivatives: 


d^' 


d^ 


^  =  i-^  +  ia(p\x  +  c.c. 


dx 


R  dJJ 


d-(p' 

dx- 


I  d'd  2ia  dd  id  da 

:t  +  — ^  ^ U  +  c.c. 


R- 


R  dx  R  dx 


(2)-(4) 


(5) 

(6) 


The  shape  function  (p  and  wavenumber  a  depend  on  the  slowly  varying  scale  x  while  the  wave 
function  x  depends  on  the  rapidly  varying  scale  x  (x  =  Rx). 


In  the  PSE  approach  the  second  derivative  temi  in  (6)  is  neglected  based  on  physical 
arguments.  Substituting  into  the  disturbance  equations  gives  a  system  of  equations  of  the  form: 
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(L„  +  i,)^  +  i,^  +  ^i,^  =  0  (7) 

OX  ox 

Here  L„is  the  compressible  version  of  the  Orr-Sommerfeld  operator  (provided  curvature  terms 
have  been  neglected),  Lj  contains  the  nonparallel  basic  state  terms,  Z^and  i^arise  due  to  the 

nonparallel  disturbance  terms.  The  system  of  equations  (7)  is  pai'abolic  and  thus  requires 
boundary  and  initial  conditions.  The  boundary  conditions  at  y  =  Oand  y  =  y^^g^^ai'e; 

u  =  v  =  w  =  T  =  0  (8) 

Notice  that  setting  Z,*  ^rid  Z3  equal  to  zero  removes  the  nonparallel  effects  and  the  problem 
(7)-(8)  reduces  to  the  linear  parallel  problem. 

There  still  remains  the  matter  of  the  ambiguity  in  x-dependence  between  ^and  x  in  the 
decomposition  (4).  This  ambiguity  is  resolved  by  imposing  any  of  the  normalization  conditions 


where  is  the  location  of  the  maximum  magnitude  of  u  (superscript  *  denotes  complex 
conjugate).  Many  other  normalizations  are  possible.  The  normalization  ensures  that  any  rapid 
changes  in  the  streamwise  direction  will  be  “absorbed”  by  the  wave  function  so  that  the  shape 
function  will  vary  slowly  in  this  direction.  This  permits  us  to  discard  the  second 

derivative  term  in  equation  (6).  Other  normalizations  may  be  used  and  will  give  slightly  different 
results.  Herbert  (1994)  is  working  on  an  “optimal  norm”  that  will  minimize  the  effect  of  the  PSE 
approximation  on  the  solution.  An  integral  normalization  may  be  used  rather  than  one  applied  at 
y  =  3’,„;,xto  avoid  the  problems  associated  with  shape  functions  developing  multiple  maxima. 

To  complete  the  problem  formulation,  initial  values  of  the  disturbance  flow  quantities  must  be 
specified  at  some  streamwise  location  (x„)  for  the  start  of  the  analysis.  If  the  analysis  begins  in  a 

region  where  the  initial  disturbance  amplitudes  are  small,  the  LST  can  be  used  to  obtain  these  initial 
conditions. 

The  nonlinear  PSE  (NPSE)  are  derived  in  a  fashion  similar  to  LPSE.  Each  disturbance 
quantity  is  transformed  spectrally  in  the  spanwise  and  temporal  directions  using 
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where 


oo  oo 


/;  =  — oo  k=- 


shapcjunction 


w'ovejmctUm 


^\n,k) 

dx 


(12) 


(13) 


Here  each  mode,  (n,k),  is  considered  to  be  the  product  of  a  shape  function  and  a  wave  function. 

Because  the  physical  disturbance  quantities  are  real, 

*  ~  ^ 


^in.k)~  ^(-,,-ky  P(n,k)  ~  P{-„-k)’  “(/i.i)  ~  ^{luk)  ~  \-n-k) 

^{n.k)  ~  ^{-n.-kp  \.k)  ~  ^(-n-k)’  \.k)  ~  ^(-n-k) 


(14) 


If  the  basic  state  is  symmetric  with  respect  to  the  ^-direction,  the  following  additional  relationships 
among  the  disturbance  quantities  hold: 


^(n.k)  ~  ^(n-kp  P{n,k)  ~  P{n-kp  %,k)  ~  ^‘{n-kp  ^(ii.k)  “  %,-k) 

'^(n,k)  ~~'^\„-k)’  ^{n.k)  -  '^{n-kp  \n,k)  ~  \n-k) 


(15) 


Substituting  (12)  into  the  disturbance  equations  gives  a  system  of  equations  of  the  form 


oo  oo 


E  E1(4,  +  L,)0  +  L2^  +  0L,^ 


(16) 


n=-ook= 


(n.k) 


The  left  hand  side  consists  of  the  lineai'  terms,  while  the  right  hand  side  consists  of  the  nonlinear 
terms.  The  portion  in  brackets  on  the  left  hand  side  contains  the  same  quantities  as  in  equation  (7) 
for  LPSE  except  the  quantities  a  and  0now  carry  the  subscripts  (n,k)  identifying  them  with  a 
particular-  mode,  and  a  and  (3  appearing  in  equation  (7)  must  be  replaced  with  n(0  o  ^nd  k(3  q 

respectively. 

The  nonlinearities  ai'e  quart ic  for  compressible  flows  where  the  effect  of  thermal  property 
fluctuations  is  considered  to  be  important.  For  incompressible  flows  this  effect  can  be  neglected 
leading  to  only  quadratic  nonlinearities.  In  fact,  it  is  reasonable  to  neglect  the  cubic  and  quartic 
nonlinearites  for  compressible  computations  at  moderate  Mach  numbers  for  cases  where  their  effect 
on  the  solution  is  small.  This  results  in  considerable  savings  since  the  higher-order  nonlinear 
terms  are  expensive  to  compute. 

Since  a  numerical  solution  to  the  system  (16)  is  desired,  a  finite  number  of  modes  {n,k)  must 
be  considered,  so  the  analysis  will  be  restricted  to  -7/  <  n<  N,  -  K  <  k  <  K .  This  is  not  a 
severe  restriction  since  we  are  interested  in  the  downstream  evolution  of  disturbances  which 
consist  of  only  a  few  modes  of  finite  amplitude  at  the  starting  location  of  the  analysis.  As  the 
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disturbances  propagate  downstream  they  may  be  amplified  and  grow  to  the  point  where  they  excite 
higher  modes  through  nonlinear  interactions.  Typically,  the  magnitude  of  the  nonlinear  terms  is 
monitored  during  the  numerical  solution  and  modes  will  be  included  in  the  analysis  as  they  become 
important. 


Harmonic  balancing  of  equations  (16)  for  a  finite  number  of  modes  leads  to  a  system  of 
equations  governing  each  (n,k}  mode  (neglecting  higher-order  nonlinearities). 


OX  ox 


\u.k) 


{,a) 


//=//,  k—k2  +ky 

li/.liv  \k,\^K  ■ 


(17) 


The  boundary  conditions  (8)  must  then  be  applied  for  all  the  modes  except  the  mean-flow 
distortion  (MFD).  This  mode  {n=k=0)  requires  a  special  boundary  condition  for  the  normal 
velocity  at  y  =  to  allow  for  changes  in  the  displacement  thickness  of  the  mean-flow  profile. 
The  boundary  condition  =  Ois  then  replaced  with 


^'(0.0)  _  Q 

dy 


(18) 


For  problems  where  there  is  basic-state  z-symmetry  we  need  only  to  solve  for  1/4  of  the  modes  {0 
<11  <  N,  0  <  k  <  K).  The  remaining  modes  will  be  required  during  calculation  of  the  nonlinear 
terms,  but  they  can  be  computed  from  the  former  set  of  modes  by  using  equations  (14)-(15).  If  no 
z-symmetiy  exists  the  computational  effort  may  still  be  halved  using  the  conditions  of  realness  of 
the  physical  disturbance  (equations  (14)).  The  nonlinear  terms  in  equation  (17)  couple  the 
governing  equations  for  all  the  modes. 

Once  the  continuous  PSE  have  been  derived,  there  are  a  number  of  discretizations  available. 
The  codes  developed  by  Bertolotti  (1990),  Stuckert  et  al.  (1993),  Stuckert  et  al.  (1994)  use 
Chebyshev  collocation  in  the  wall-normal  direction  and  backward-Euler  finite  differences  for  the 
streamwise  direction.  The  backward-Euler  finite  difference  helps  damp  transients  that  occur  as  a 
result  of  approximations  used  to  obtain  the  initial  conditions.  Bertolotti  (1990)  uses  backward- 
Euler  discretization  near  the  initial  location  and  slowly  switches  to  Crank-Nicolson  further 
downstream.  Chang  et  al.  (1991)  use  either  a  fourth-order  central-difference  or  compact  two-point 
scheme  for  derivatives  in  the  wall-normal  direction  and  second-order  backward  differences  (multi- 
step)  in  the  streamwise  direction.  Haynes  &  Reed  (1996)  use  fourth-order  central  differences  for 
derivatives  in  the  wall-normal  direction  and  backwai'd-Euler  discretization  in  the  streamwise 
direction. 


Reed:  Hypersonics  (AFOSR) 


page  14 


Although  the  removal  of  pai't  of  the  second-derivative  terms  in  the  PSE  formulation  apparently 
results  in  a  system  of  parabolic  differential  equations,  there  is  still  some  elliptic  behavior  associated 
with  the  upstream  propagation  of  acoustic  disturbances.  This  situation  is  analogous  to  the 
parabolizing  procedures  used  to  develop  the  PNS  equations  (Anderson  et  al.  1984).  For 
incompressible  formulations  the  streamwise  derivative  of  the  shape  function  for  pressure  may  be 
dropped,  or  a  lai'ge  enough  stepsize  can  be  used  to  “step  over”  the  elliptic  region  (Malik  &  Li 
1993).  The  two-dimensional  incompressible  formulation  of  Bertolotti  (1990)  avoids  this  difficulty 
by  using  a  streamfunction  formulation.  For  compressible  PSE  the  resulting  equations  ai'e 
hyperbolic  in  the  supersonic  region  and  elliptic  in  the  subsonic  region.  Here  the  Vigneron 
technique  (Vigneron  et  al.  1978)  can  be  used  to  suppress  the  upstream  wave  propagation  in  the 
subsonic  portion  of  the  boundary  layer. 

Wc  have  formulated  and  coded  the  PSE  formulation  for  a  surface  (including  curvature).  Then 
wc  verified  and  validated  our  code.  The  basis  of  validation,  or  confirming  that  the  equations  used 
to  model  the  physical  situation  are  appropriate,  is  assumed  to  be  a  successful  comparison  with  the 
few  careful,  archival  experiments  available  in  the  literature. 

Over  the  next  year,  equilibrium  effects  will  be  included  and  the  flow  over  the  circulai'  cone 
investigated  using  the  PSE  analysis. 

Verification 

We  consider  verification  to  mean  “confirming  the  accuracy  and  correctness  of  the  code”. 
There  are  mainly  three  sources  of  error  in  the  abstraction  of  continuous  PDE's  to  a  set  of  discrete 
algebraic  equations;  (1)  discretization  errors,  (2)  programming  errors  (bugs),  and  (3)  computer 
round-off  errors.  Of  the  above  three,  only  programming  errors  can  be  completely  eliminated.  The 
objective  of  code  verification  is  then  to  completely  eliminate  programming  errors  and  confirm  that 
the  accuracy  of  the  discretization  used  in  solving  the  continuous  problem  lies  within  some 
acceptable  tolerance.  Aside  from  specifying  single  or  double  precision,  the  code  developer  has  little 
control  over  the  computer  round-off  errors,  but  this  is  usually  several  orders  of  magnitude  smaller 
than  the  discretization  error  and  far  less  than  the  desired  accuracy  of  the  solution. 

In  this  section  we  address  programming  and  discretization  errors.  Many  methods  ai'e 
discussed  in  the  literature  for  code  verification  using  grid  refinement,  comparison  with  simplified 
analytical  cases,  etc.  For  recent  discussions  see  Roache  (1997)  and  Oberkampf  et  al.  (1995). 
Specific  suggestions  for  testing  a  CFD  code  for  the  study  of  transition  include  (a)  grid-refinement 
studies,  (b)  solving  test  problems  for  which  the  solution  is  known,  (c)  changing  the  “far-field” 
boundary  locations  systematically  and  re-solving,  (d)  comparing  linear  growth  rates,  neutral 
points,  and  eigenfunctions  with  linear  stability  theory,  (e)  running  the  unsteady  code  with  time- 
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independent  boundary  conditions  to  ensure  that  the  calculations  remain  steady,  and  (f)  running 
geometrically  unsymmetric  codes  with  symmetric  conditions. 

In  addition  to  the  usual  code  verification  techniques,  there  is  a  general  method  which  we  used 
to  verify  the  discretizations  and  locate  programming  errors  by  comparison  with  “manufactured” 
analytical  (Steinberg  &  Roache  1985).  This  method  is  general  in  that  it  can  be  applied  to  any 
system  of  equations.  Although  it  is  an  extremely  powerful  tool,  this  method  has  received  little 
acclaimation  in  the  literature. 

For  clarity  the  technique  is  demonstrated  on  the  Poisson  equation. 

,  d~u  d^u  .  , 

LM  =  ^  +  ^  =  F(x,y)  (19) 

dx  dy 

To  solve  this  problem  we  discretize  the  operator  L  using  some  appropriate  approximation  (finite 
differences,  spectral,  etc.).  In  general,  the  exact  solution  is  not  available.  Therefore,  for 
verification  purposes,  wt  force  the  solution  to  (19)  to  be  some  combination  of  analytical  functions 

with  nontrivial  derivatives.  For  example,  we  might  consider  the  system  g  =  Lv  =  sin(2x), 

which  has  an  analytical  solution  v  =  sin(2x).  The  exact  solution  can  then  be  compared  with 
the  computed  solution.  Of  course,  manufactured  solutions  should  be  chosen  with  topological 
qualities  similai'  to  those  anticipated  for  the  solution  to  the  “real”  problem  (e.g.  gradients  close  to 
the  wall).  Proper  choice  for  the  manufactured  solutions  also  allows  the  discretization  of  the 
boundary  conditions  to  be  verified.  For  large  systems  of  equations  a  symbol  manipulator  is 
recommended  for  computing  g.  If  a  bug  occurs,  zeroing  the  coefficients  of  some  terms  in  equation 
(16)  can  help  to  isolate  the  bug. 

We  recommend  this  method  in  addition  to  the  other  specific  suggestions  mentioned  above  for 
code  verification.  Code  validation  is  discussed  in  the  next  Section. 

Validation 

To  date  the  PSE  have  been  applied  to  a  variety  of  2-  and  3-D  flow  situations  and  are  generally 
regai'ded  as  appropriate  for  convectively  unstable  flows  (Haynes  &  Reed  1996;  Schrauf  et  al. 
1995;  Stuckert  et  al.  1994;  Wang  et  al.  1994;  Stuckert  et  al.  1993;  Herbert  &  Lin  1993;  Herbert 
1994;  Malik  &  Li  1993;  Bertolotti  et  al.  1992;  Chang  et  al.  1991 ;  Bertolotti  1990). 

Bertolotti  et  al.  (1992)  verified  the  PSE  approach  for  T-S  (two-dimensional)  disturbances  in  a 
Blasius  boundary  layer  by  comparison  with  the  DNS  results  of  Bayliss  et  al.  (1986).  For  this 
nonlinear  case  they  compared  not  only  the  growth  rates,  but  also  the  mode  shapes  of  the  harmonics 
and  found  excellent  agreement.  Three-dimensional  nonlinear  PSE  stability  results  were  compared 
with  the  experimental  results  of  Kachanov  &  Levchenko  (1984),  but  only  qualitative  agreement 
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was  achieved.  The  differences  are  attributed  to  virtual  leading-edge  and  slight  pressure-gradient 
effects  in  the  experiment.  Comparison  of  the  same  PSE  results  with  DNS  results  of  Easel  et  al. 
(1990)  and  Crouch  (1988)  show  better  agreement.  The  experimental  results  of  Cornelius  (1985) 
for  K-type  transition  are  compared  with  PSE  results  and  qualitative  agreement  is  found. 

All  of  the  above  results  are  on  a  flat  plate.  The  only  available  experiments  against  which  to 
validate  our  3-D  cuiwed-suiface  code  are  those  of  Reibert  et  al.  (1996)  on  an  NLF(2)-0415  swept 
airfoil.  The  nonlinear  PSE  results  include  curvature  effects.  It  is  clear'  that  the  linear  theories  fail  to 
accurately  predict  the  transitional  flow  for  this  situation.  Studying  a  comparison  of  the 
experimental  and  computational  disturbance  mode  shapes,  we  demonsti'ated  that  the  nonlinear  PSE 
does  an  excellent  job  of  capturing  the  detailed  profiles.  See  Haynes  &  Reed  (1996). 

Global  Eigenvalue  Solver 

To  generate  initial  conditions  for  the  nonlinear  PSE  code,  we  have  also  developed  a  global 
eigenvalue  solver  to  predict  which  disturbances  ai'e  most  likely  to  break  down  to  turbulence.  We 
apply  this  code  to  the  results  of  the  basic-state  code,  and  then  apply  the  PSE  to  predict  transition 
location. 
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